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I \ Abstract. We study probability distributions of eigenvalues of Hermitian and non- 

^ ■ Hermitian Euclidean random matrices that are typically encountered in the problems 

of wave propagation in random media. 
^ ■ 

I ! PACS numbers: 42.25.Dd, 02.50.Cw, 05.40.-a 

2 ■ 1. Introduction 

^ \ Random matrix theory is a powerful tool of statistical physics [Ij with important 

^ i applications in the field of quantum and wave transport in random media [21 [3l [4]. 

f-^ \ A special class of random matrices are Euclidean random matrices with elements Fij 

^^ ' defined with the help of some function /(r^, Vj): F^j = /(r^, r^). Here r^ (i = 1, . . . A^) 

are randomly chosen points in the Euclidean space [51 [6] . Euclidean random matrices 
O \ appear in various physical contexts and were previously considered to interpret the 

'boson peak' in supercooled liquids |7j or to study slow relaxation in glasses and scalar 
phonon localization [8], to cite a few recent examples. The purpose of this paper is to 
k>^ \ study eigenvalue distributions of certain large Euclidean random matrices that appear 

^ \ in problems of wave propagation in random media. Because in the simplest case of 

scalar waves the propagation is described by a scalar wave equation, the function / that 
will be of interest to us is the Green's function G(r^, r^) of the Helmholtz equation 

Ajr 

(V^ + kl + iv) G(r„ r,) = --5(r, - r,), (1) 

/Co 

where 77 is a positive infinitesimal. It is easy to check that G(r^,rj) = exp(i/co|r^ — 

r^-|)/A;o|r^-r^-|lj|. 

Statistical properties of the ensemble of matrices G with elements Gij = (1 — 

Sij)G{ri^ Yj) for A^ ::> 1 are of primary importance in the context of Anderson localization 

of electromagnetic [9l [10, [H] and matter [12] waves. The same matrices appear in the 

studies of collective spontaneous emission in dense atomic systems [131 El tlSl [IS [Ej . 

i We restrict ourselves to three-dimensional space in this paper. 
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The interplay between Anderson localization and Dicke superradiance can also be 
described by this ensemble of matrices flS] and properties of their eigenvalues are 
important for understanding of random lasers p!9l [20] and dynamic instabilities in 
nonlinear random media |21|. Meanwhile, the matrix G is non-Hermitian, its eigenvalues 
are complex and their probability distribution is difficult to access. This is why in 
several works dealing with superradiance [HI [HI [TTl [18] the imaginary part of (7, a 
matrix with elements sin(/co|r^ — rj|)//co|r^ — rj|, was considered. This real symmetric 
matrix is much easier to study and in many situations it still contains some of the 
important aspects of the full problem. Similarly, the real part of (7, a matrix with 
elements cos(/co|r^ — rj|)//co|r^ — rj|, is relevant for understanding the collective Lamb 
shifts in dense atomic systems [l5l [IT] . 

Despite the importance of the three matrices G^ S = ImG and C = ReG introduced 
above little is known about statistical properties of their eigenvalues. In the general 
case, the eigenvalue distribution of G was studied only numerically P [IHl [H]. Some 
analytic results are available in the limit of high density of points r^ inside a sphere: 
p = N/V -^ oo [m [151 [151 [II] 5 when the summation in the eigenvalue equation 
V . Gijipj = Xipi can be replaced by integration. The purpose of this paper is to 
partially fill this gap by considering eigenvalue distributions of the three matrices above 
at finite densities p, with the distances between neighboring points r^ that are larger 
than, comparable, or smaller than the wavelength Aq = 27r//co. This situation is of 
particular importance in the context of wave propagation in random media because in 
order to observe phenomena due to scattering of waves on the heterogeneities of the 
medium, the density of scattering centers (or scatterers) should be neither too low (in 
this case the scattering is negligible), nor too high (in this case the medium responds 
as an effective homogeneous medium). One of the possible experimental realizations 
of a strongly scattering system is a cloud of cold atoms in which propagation of quasi- 
resonant light (wavelength Aq) is studied [§|. Nowadays such clouds are routinely created 
at densities pAg <C 1, allowing observation of interesting phenomena due to the multiple 
scattering of light [22l [23]. This justifies the importance of properly understanding the 
low-density regime. However, the most interesting phenomena for waves in an ensemble 
of point-like scattering centers are known to take place at densities pAg ^ 1, when 
interference effects become important, eventually leading to Anderson localization (see, 
e.g., [24l [25l [26] and references therein). Our results may be useful for understanding 
Anderson localization and its interplay with other collective phenomena (such as Dicke 
superradiance) fiE\ . 



Light is a vector wave but here we restrict ourselves to a scalar approximation. 
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2. Summary of main results 

Before presenting the details of calculations, let us list our main results: 

• A general framework is developed to deal with Hermitian Euclidean matrices 
(section [3]) . We show how the theory of asymptotically free random variables can 
be applied in this context. 

• The approach developed in section [3] is applied to study the probability distribution 
^(A) of real eigenvalues A of the real symmetric random matrix S corresponding to 
A^ points in a box of side L (section [4]) . We show that when /3 = 2.8N/{koLy < 1, 
p{X) is given by the famous Marchenko-Pastur with /3 = varA as the only parameter. 
For /3 > 1, the Marchenko-Pastur law does not apply anymore. 

• The probability distribution p{X) of real eigenvalues A of the real symmetric random 
matrix C is studied (section Ej) . We show that p{X) depends on two parameters: 
/3 and the number of points per wavelength cube pAg. Analytic results are in 
agreement with numerical simulations for pAg ^ 30 and any /3. In the low-density 
limit pAg <C 1, p{X) exhibits a transition from the Wigner semi-circle law for /3 <C 1 
to the Cauchy distribution for /3 ;» 1. 

• As the first example of non-Hermitian Euclidean matrices, in section [6] we study the 
complex symmetric matrix X = C + i{S^ — I), where two different and independent 
sets of points {r^} and {r^} are used to define the matrices C and S\ For /3 < 1, 
the probability distribution of complex eigenvalues A of X is obtained by combining 
the results for S and C obtained in section [4] and section [5l respectively, in the 
framework of the theory of free random variables. The domain of existence of 
eigenvalues of X undergoes a transformation from a circular to a triangular shape 
as (3 increases from to 1. For /3 ;» 1, numerical simulations show that the support 
of the distribution on the complex plane takes an 'inverted T' shape. 

• The non-Hermitian matrix G = C + i{S — I) with elements given by the Green's 
function of the Helmholtz equation ([T]) is studied in section [71 by means of extensive 
numerical simulations. We find that at low density pAg ^ 30 and for /3 <C 1 
the domain of existence of eigenvalues of G on the complex plane coincides with 
that of X and is given by a circle of radius ^/2^ centered at (0, |/3). At larger /3, 
the domain remains approximately a circle with the same center but a larger radius 



R ^ \/ 2/3 + (|/5)^. When the density pAg reaches a critical values of approximately 
30, a 'hole' opens in the eigenvalue distribution that otherwise still keeps its circular 
shape. 

The numerically evaluated marginal distributions of real and imaginary parts of the 
eigenvalues A of the matrix G roughly follow the laws obtained for the eigenvalues 
of the matrices C (for pAg ^ 30) and S (for |/3 < 1), respectively. For i/3 > 1, the 
distribution of F = ImA + 1 approaches the 1/F law. 

The mean minimum value of ImA is approximately given by (min(ImA)) o^ —1 + 
2.3/(A X pAg)^/^ for pX^ < 10 and decays faster at higher densities. The mean 
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maximum value of ImA is roughly (max(ImA)) c^ ^/3 + R. 

The above mathematical results have important applications in a number of 
physical problems of contemporary interest, as discussed in section El In particular, 
they provide an additional insight into the cooperative spontaneous emission of large 
atomic clouds (section 18. ip . Anderson localization (section 18.21) . and random lasing 
(section 18.31) . 

3. General framework 

Consider a singly-connected three-dimensional region of space V. Let {'0^(r)} be an 
orthonormal basis in V, such that 

dh ^^{r)rn{r) = S,r^n• (2) 

V 

We will now show that an arbitrary N x N Euclidean random matrix F with elements 
F,, = /(r„r,), iJ = l,...N, (3) 

where / is a sufficiently well-behaved function of r^, Yj G V^ can be represented as 

F = HTHl (4) 

Here H is di N x M matrix with elements 

Him = J—i^m{ri). (5) 

We use V to denote the considered three-dimensional region of space as well as its 
volume, T is a M X M matrix to be defined below, and the dagger 'f denotes Hermitian 
conjugation. The size M of the matrix T can be arbitrary and, in fact, M will be infinite 
for the majority of functions /(r^, r^). 

To establish (|4j), we write the zj'th element of the matrix F explicitly as 

F^, = ^J2^mnM^^)rnirJ), (6) 

where we used ([5]) and the definition of matrix multiplication. Multiplying this equation 
by '0^/(r^)'0^/(rj), integrating over r^ and r^, and using the orthogonality of the basis 
functions '0m(r), we readily obtain 

Tmn = yjdh,ldhj /(r,,r,-)C(ri)^n(r,-). (7) 

V V 

It is easy to check that with the elements T^^ of T defined by ([7]), (|4]) is indeed obeyed. 

When the points {r^} are chosen inside V randomly, F and H become random 

matrices, whereas T is always a non-random matrix independent of {r^} and determined 

uniquely by the function /, the region V^ and the choice of the orthonormal basis 
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{ipm{r)}. We will limit our consideration to the case when the spatial integral of any 
basis function iprn{^) that contributes to (|6]) vanishes [|||: 

dhiljm{r) = 0. (8) 

V 

The elements Hi^ of H are then independent random variables having zero means and 
variances equal to 1/A^: 

(H.^) =\^ j d^r, \l^^m{vi) = 0, (9) 

V 



(HrmH;^) = ^Jd'r,Jd'rj ^^^(r,)C(ri) 



V V 

= {H,m){H;^)=0, % + 3. (10) 

^H^miil.) = l^jdh, ^i^m{r^)rn{r^) = ^. (H) 

V 

The representation (|4]) is very useful because it can be dealt with using the powerful 
mathematical arsenal of the so-called free random variable theory [271 EH EH] . Without 
going into details, we remind the reader that for random matrices, the notion of 
asymptotic freeness [27] is equivalent to the notion of statistical independence that we 
are familiar with for random variables. Three fundamental objects of the free random 
variable theory, defined for any Hermitian matrix F, will be useful for us in this paper: 
the usual Green's function 

the Blue function B{z) equal to the functional inverse of Q{z): 

B[g{z)] = z, (13) 

and the S'-transform of the probability distribution of eigenvalues defined through an 
auxiliary function x(^)- 

S{z) = ^^^xiz), (14a) 

Z 



1 



x{z)' 



1 



x(^). 



1 = z. (146) 



If two Hermitian random matrices A and B are asymptotically free, the Blue function 
B^{z) of their sum (7 = A + S is equal to the sum of individual Blue functions B^{z) 
and B^{z)^ minus 1/ z. The S'-transform of the matrix product C = AB can be found 
by multiplying the individual S'-transforms of A and B. Once the Blue function or the 
S^-transform corresponding to the random matrix C are found, its Green's function Q{z) 

II This restricts the class of functions /{vi^Vj) to which our analysis applies but will be sufficient for 
us here. 
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6 



can be calculated either from ( IT31) or from ( I14al) and ( llAbh . The probability density of 
the eigenvalues A of (7 is then determined in the usual way: 

1 



p{X) = — — lim Im^(A + ie). 



(15) 



The functions G{z)^ B[z) and S{z) all contain the same full information about the 
statistical distribution of eigenvalues A as p{\). The Green's function can be represented 
as a series with coefficients in front of consecutive powers of 1/ z equal to statistical 
moments of A: Q{z) = Yl'^=o{^^) / ^^^^ - We have, therefore, 

1 rf^+i^(z) 
(n + l)!rf(l/z)^+i 

where z is assumed real. Using this equation and ( TT31) we readily derive an expression 
for (A^) in terms of B{z): 



(A") 



(16) 



(A^ 



1 



(n + 1)! 



' B'{z)dz 



B'{z) 



(17) 



where B\z) = dB{z)/dz. If we introduce the /^-transform TZ{z) = B{z) — \j z [28j, 
the average eigenvalue and the variance become (A) = 7^(0) and varA = ((A — (A))^) = 
7^' (z) 1 2^0 5 respectively. 

For matrices F of the form dH) , the free random variable theory provides a number 
of mathematical theorems that we will exploit in this paper. In particular, one shows 

that 

1 



^^^^^"z + m/tv^^Im^ 



(18) 



if T is a Hermit ian nonnegative random matrix independent of H and the limits A^, 
M ^ oc are taken at a constant M/N . Using ( IT81) . we derive a relation between the 
Blue function of F and the Green's function of T: 



-^<^)4{-# 



ri 


fA 


' 


-Qf 




-1 


[z 


\zj 





(19) 



A particular case that we will consider in the remainder of this paper is when the 
region ]/ is a square box of side L [see figure [l](a)]. A convenient set of basis functions 
is then given by 'plane waves' 

1 



'^mW 



:^iqm-r 



(20) 



where q^ = {<?m.., ^m^, ^mj, ^m. = rn^^q with m^ = ±1, ±2, . . . (and similarly for q^y 
and qrriz)) ^^d Ag = 2ti/L. Equation ([71) is then simply a double Fourier transform of 
the function /(r^, r^) in the box and the representation (|4j) stems from the Fourier series 
expansion of /(r^,rj), without the harmonics corresponding to q^ = 0. 
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N points 





27ra ^^ 27ra^ 
L ^^ L 



Figure 1. (a) We consider TV points randomly distributed in a three-dimensional 
cube of side L. (b) Regions in the Fourier space. For the sine random matrix, only 
the region 2 contributes to the matrix T. In contrast, for the cose random matrix, T 
has contributions from the regions 1 and 3 but not from the region 2. 



4. Eigenvalue distribution of the sine matrix 

We start by considering the real symmetric N x N Euclidean matrix F = S with 
elements defined through the cardinal sine (sine) function: 

sin(/co|r^ -r^-|) 



Sij — f\Vi Yj) 



fcol^i ~ ^j\ 



(21) 



Here /cq is a constant and the vectors r^ define positions of A^ randomly chosen points 
inside a three-dimensional cube of side L. 

The first important property of the matrix S is the posit iveness of its eigenvalues: 
A > 0. Indeed, the Fourier transform of the function /(Ar) in f l2T]) is positive and hence 
/(Ar) is a function of positive type. An Euclidean matrix defined through a function 
of positive type is positive definite and hence has only positive eigenvalues. The matrix 
T corresponding to S can be found from ([7]): 



T 

■J- rr 



-v^h'-l 



d^ra 



sin(A;o|ri -ra]) 
kQ\^i -r2| 



-iqmri+iq„r2 



(22) 



V V 

Unfortunately, it is impossible to calculate this double integral exactly in a box. 
However, introducing new variables of integration R = |(ri +r2) and Ar = r2 — ri and 
limiting the integration over Ar to the region |Ar| < L/2a^ with a ^ 1 a numerical 
constant to be fixed later, we obtain an approximate result 



-J- m.n. — 



= 5r7 



A^ 



(j3j^ g-i(qm-qn)R 



d^Ar 



sin(/coAr) 
koAr 



,i(qm+qn)Ar/2 



V 



27r^N L 



smc 



{qn 



|Ar|<L/2a 



smc 



{Qn 



^4' 



(23) 



' koq^V 2a7i 

This expression is still too involved to be useful. In order to simplify it, we note that 
the second sine function in (1251) is always smaller than 2a/kQL (because Qm = |qm| > 
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and ko > 0) and hence can be dropped in the hmit of large k^L ^ 1 considered in this 
paper. Furthermore, because the first sine function in f l23l ) is peaked around Qm = ko^ 
we replace it by a boxcar function n[(g^ — ko)L/2a7r]^ where n(x) = 1 for |x| < | and 
n(x) = otherwise. The coefficient in front of (g^ — ko) in the argument of 11 is chosen 
to ensure that the integral of the latter over g^ from to oc is equal to the same integral 
of the sine function. We then obtain 



27r^N L ^ 

T ^ n 

"^" - e,V 2a7r 



{qm - ko)- 



(24) 



27Ta 

which is different from zero only for q^'s inside a spherical shell of radius /cq and thickness 
27Ta/L [i.e. in the region 2 of figure [T](b)]. In addition, for all q^'s inside the shell the 
value of Tmn is the same and equal to N/M with M = a^koLY/n ^ 1 the number of 
q^'s inside the shell. Equation (|4]) then yields 

S = —HH^ (25) 

which is equivalent to (|lj) with a M x M matrix T = {N/M)l^ where I is the identity 
matrix. We then readily find gf{z) = {1/M)Tt[z - (N/M)!]-^ = (^ - N/M)'^ and 
from ([I9j): B^{z) = (1 - /3^)"^ + l/z with /3 = N/M. This is the Blue function of the 
famous Marchenko-Pastur law [28l [30] : 

p{X) = (l - i)^(A) + ^Z" ' A„.„K(A„„ ^AK^ pgj 

where Amin,max = (1 T V^)^ and x+ = max(x, 0). The distribution of eigenvalues of the 
matrix f l2Tl ) is therefore parameterized by a single parameter (5 equal to the variance of 
this distribution, as it is easy to check from f l26l) : var(A) = /3. 

Although we derived f l26l) using the machinery of free random variables applied to 
Euclidean matrices as discussed in section [3l it represents a somewhat trivial example 
of application of this technique because the matrix T in (|4j) turns out to be proportional 
to the identity matrix I. Equation f l26l) was first derived long before the theory of 
asymptotically free random variables was introduced [30j. It can be established using 
various approaches, such as, e.g., the diagrammatic technique [31j. However, to our 
knowledge, the fact that this distribution describes eigenvalues of the Euclidean matrix 
S was never noticed before. The advantage of using the free random variable theory to 
study Euclidean random matrices is that f l26l ) now appears as a special (and apparently 
the most trivial) case of a wide class of distributions describing matrices of the form (|4j) . 

Note that despite the fact that our derivation of f l26l ) was based on several 
approximations, the average value of A, (A) = 1, following from this equation is exact. 
The second moment of A can also be found directly from ( l2TI) . For /cqI/ 3> 1 and in the 
limit A^ ^ oc we find: 

where the numerical constant a is given by 

a = \ f d^ui [ dW, ^- ^^2.8, (28) 

2 J J U1-U2M 
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1.0 



numerical 

Marchenko-Pastur law (26) 
N = 10"^ 

(3 = 0.33 




12 3 4 

Eigenvalue A 



5 10 15 

Eigenvalue A 



Figure 2. Probability density of eigenvalues of a square N x N Euclidean matrix S 
with elements Sij = sin(/co|r^ — rj|)//co|r^ — rj|, where the N points r^ are randomly 
chosen inside a 3D cube of side L. Numerical results (blue solid lines) obtained for 
A^ = 10^ after averaging over 10 realizations are compared to the Marchenko-Pastur 
law (j26j) (red dashed hues) with /3 = 2.8N/{koL)'^ for several densities p of points 
(Ao = 27r/A:o). 



with the integrations running over the volume of a cube of unit side. By requiring that 
the second moment 1 + /3 of the distribution f l26l) coincides with f l27j) we can now fix the 
value of a that remained arbitrary until now. We obtain a = n/a c::^ 1.12 and 

2.8A^ 



/3 



(29) 



In figure [2] we present a comparison of f l26l ) with the results of direct numerical 
simulations. The latter amount to generate A^ random points r^ inside a three- 
dimensional cube, to use these points to define a random N x N matrix S according 
to ( 1211) . and to diagonalize S using the standard software package LAPACK [|32j. The 
procedure is repeated several times and a histogram of all eigenvalues A is created. This 
histogram approximates the eigenvalue distribution p{X). As we see from figure O the 
agreement between numerical results and the Marchenko-Pastur law f l26l ) is good for 
/3 < 1 but f l26l) fails to describe p{X) when (3 becomes larger than unity. The reason 
for this is easy to understand if we go back to fl22l). f l23l ) and f l24l ). Indeed, when we 
approximate the result of integration in f l22l ) by f l24l ) , we reduce the infinite-size matrix 
T to a matrix of finite size M x M. By definition, the rank of the latter matrix is inferior 
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or equal to M. The rank of S = HTW cannot be larger than the rank of T and hence 
is also bounded by M from above when we use f l24l ). When /3 > 1, implying M < N ^ 
the representation (|4]) only gives us access to M of A^ eigenvalues of S', which is not 
sufficient to reconstruct the probability density p{\). In order to access the regime of 
/3 > 1 one needs to find a better approximation to f l22l ) than f l24l ). 

Note that the eigenvalue distribution of the matrix S has been studied numerically 
by Akkermans et al. in the context of light propagation in atomic gases (see figure 
1 of [18]) without proposing any analytical approximation to it. The parameter 
(5 ^ N/{koLy has been introduced in that work as a ratio of the number of atoms A^ to 
the number of transverse optical modes A^^ oc (/cqI/)^. The same parameter appeared in 
[HI [El [ini [17] as a superradiant decay rate in a cold atomic gas. Hence the results of 
this section complement and extend the works [HJ [151 [iHl [13 [H] . 

5. Eigenvalue distribution of cose matrix 

Let us now consider an Euclidean random matrix with elements defined using the 
cardinal cosine (cose) function: 

Q, = /(r. - r,) = (1 - S^, rf;^'^-f\ (30) 

^0 I A i — Tj I 

The prefactor 1 — 5ij allows us to deal with the divergence of the function cos(x)/x for 
X ^ 0. However, in the beginning of our analysis we will ignore this prefactor and will 
use Cij = cos(/co|r^ — rj|)//co|r^ — Tj\ for all i^j (the diagonal elements of C are thus 
infinite). Proceeding as in the previous section, we find 

AttN 1 
koV ql-kl 

under the same approximations as in ( l23l) (i.e., we extended integration over Ar to the 
whole space). The matrix T defined by (IHTI) has infinite size. 

The divergence of f lSB for g^ -^ ko can be traced back to the neglect of the 
finiteness of the volume V when extending integration over Ar to the whole space. 
Taking into account the fact that Ar cannot exceed a maximum value of the order of L, 
the divergence is regularized and the resulting T^^ changes sign rapidly but continuously 
in a strip of width ^ 1/L around g^ = ko. In the following, we will neglect the 
contribution of Qm's inside the spherical shell corresponding to this strip because (i) the 
shell has small thickness in the limit of koL ^ 1 that we are interested in and (ii) Qm's 
situated symmetrically with respect to the surface g^ = ko yield contributions of roughly 
equal magnitudes but opposite signs which approximately cancel. More precisely, we 
will exclude a shell of thickness 2Tia' / L around q^ = ^o and will use f lHTl ) outside this 
shell [see figure [l](b)]. The numerical constant a' r^ 1 will be fixed later. The matrix C 
therefore takes the form: 

C = -H^^^f^^m^^^^ + i/(3)r(3)i/(3)t. (32) 



-^mn — 7^ T7 ^9 u2^^ W / 
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Here the first term describes tlie contribution of qm < ko — a^n/L with Trnn = 
4:7rN/[koV{kl — q'^)]Smn- The matrix T^^"^ is a diagonal square matrix of size Mi o^ 
(47r/3)(fco — QfV/L)^(L/27r)^ obtained by dividing the volume of a sphere of radius 
ko — a'n / L [region 1 corresponding to the inner sphere in figure [11(b)] by the volume 
(27r/L)^ associated with a single mode. The second term in f l32l ) corresponds to 
9m > ^0 + QfV/L [region 3 in figure [11(b)] . The matrix T^^^ is, again, diagonal, with 
elements Trnn = 4:7rN/[koV{q'^ — kl)]5jnn but, in contrast to T^-^^ has infinite size. We will 
treat this matrix as a finite-size matrix of size Ms — (47r/3)[g^^^ — (A;o+aV/L)^](L/27r)^, 
corresponding to taking into account only g^ < ^max- The limit of g^ax ^ co will be 
taken at the end. The minus sign in front of the first term in f l32l ) was introduced to 
work with a positive-definite matrix T^^\ 

The Green's function of the matrix T^-^^ is 

1 ^ 1 



QfiD (z) = 



1 



E 



Ml ^-^ z 

qraKko—a'-K/L 



koV kl-ql^ 



Mi{p\l) 



fegi 



d«: K 



1 



PA3. 



(33) 



27r2 1-/^2 

where the wavelength Aq = 2ti /k^ and p\^ is the average number of points r^ per 
wavelength cube. The last line of this equation was obtained in the limit of k^L ^ 1 
by approximately replacing the summation over a set of discrete wavevectors q^ by an 
integration over k = Qm/^o- The integral in f l33l) can be evaluated yielding 

6^^(i)(^) = ^^^<! f 1-^1 .^^fl-"" 



ttMiZ 



koL 



3pXl 



koL 



ML 

27r^z 



1 arctan 



-1 a 7T 

koL 



27T^Z 



(34) 



1 



A similar calculation can be performed for the Green's function of T^^' except that 
the integration in (l33l) extends from 1 + a'li/koL to ^maxi Mi is replaced by M3, and 
1 — K^ in the integrand of flHHl) — by k^ — 1. We find 



Of (3) (Z) 



2N 

TvMsZ 

1 + 



Hv-n 



a'n 



kf)L 



1 

+ - 

z 



Mi 

27t'^z 



a'n 
koL 

max ~r J- ~r 

arctan 



2w^ 
3pXl 



Kj^ 



0L1\ 



koL 



1 

+ - 

z 



1 + ^^ 



Mi 

'27T^Z 
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arctan 



Kjt} 



21:^ z 



- 1 



(35) 



Because the two terms in f l32l ) correspond to contributions of different parts of q- 
space, they are asymptoticaUy free and hence the Blue function of their sum (i.e. of 
the matrix C) can be found as a sum of their respective Blue functions. The Blue 
functions of H^^)f^^)H^^)^ and H^^)f^^)H^^)^ are found using ([I9D, whereas the Blue 



function of 
sum C = — 

l/z: 



_^(i)7^(i)^(i)t ig equal to -B^wfWHmi-z). The Blue function of the 
fj{i)f{i)fy{i)] ^ fj{3)f{3)-fj{3)^ -g ^ g^^ Qf individual Blue functions, minus 

-^c(^) = ~ -SH(i)f(i)^(i)t(~^) + Bj^(3)f(3)^(3)f{z) — l/z 



1 + ^ 



TT 



arctan 



arctan 






+ 



TV 



1 + 



pAg 



Ml 
1- 



pAg 

27r2- 



arctan 



pAg 

87r/3' 



-1 



pAg, 

27r2' 



Kv| 



-1- ^z 



(36) 



where /3^ = ^N/a\k^Lf. 

The final step consists in taking the limit ^^ax -^ co. We now recall that up to now 
we ignored the fact that the matrix C had zero diagonal elements Cu = 0. Instead, we 
considered a matrix with infinitely large diagonal elements. Such a matrix naturally has 
infinite eigenvalues and to go back to the case of Cu = we have to shift the eigenvalues 
to the left. To determine the exact shift, we compute the average of A from ( l36l) using 
( ITTj) and subtract it from f l36l) because we know that for the matrix C defined by f l30l) . 
(A) = {1/N){TtC) = exactly Ij. We then compute the second moment (A^) and 
require that its value in the limit of pXl//3' oc l/k^L ^ is equal to /3 defined by f l29l ). 
This fixes (5' = (7r^/4)/3 corresponding to a^ ^ 0.45. The final expression for the Blue 
function of C is 



B, 



c 



Z TV pXq TT 



-1 



arctan 



1 + 



pAg 

2^3/3 



pA| 

27r2' 

1- 



_1 _ P^o ~ 

^ 27r2^ 



— arctan 



pAg 

2^3/3 



-1- ^Z 



TV 
2 



(37) 



The Green's function Gc{z) can be found from this equation by solving B^[Q{z)] = z 
which, for the general case, we do numerically. 

Let us consider the low-density limit of f l37l) . pAg <C 1. For large box size L ^ l/ko 
the arguments of arctan functions in ( 1371) are close to — i. They can be thus expanded 

^ Subtracting a constant from the Blue function B{z) results in shifting the eigenvalue distribution 

p{\). 
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Figure 3. Probability density of eigenvalues of a square N x N Euclidean matrix C 
with elements Cij = (1 — Sij) cos{ko\ri — rj|)//co|r^ — r^l, where the N points r^ are 
randomly chosen inside a 3D cube of side L. The left panel corresponds to the low- 
density limit and is obtained using (|38|) with /3 = 0.1, 0.5 and 5. The distributions are 
symmetric and vanish for |A| > A* with A^,. given by (|39|) . The right panel illustrates 
our equation (|4Q|) obtained in the high-density limit for two densities pAg = 20 and 50. 
For pAg > 30.3905 the distribution develops a gap in between Ai and A2 given by dH]) 
and (|^2]) . respectively. 



in series in the vicinity of this point. In the resulting expression we take the limits of 
pAo -^ and pAg//? ^ l/koL ^ to obtain 

1 1 l-f/3z 
In ^ 

z 



Bci^) 



llni 



pXl < 1. 



(38) 



This expression has two important limits. For /3 <C 1 we find B^{z) = f3z + 1/z 
which is the Blue function of the Wigner semi-circle law p{X) = ■\/4:/3 — X^ /2ti[5. In the 
opposite limit of /3 :» 1 we have B^{z) = — i + 1/z, which corresponds to the Cauchy 
distribution p{\) = l/[7r(l + A^)]. Equation ( l38l) therefore describes a transition from 
the Wigner semi-circle law at /3 <C 1 to the Cauchy distribution at /3 -^ oc. The 
eigenvalue distribution following from f l38l ) is always symmetric with respect to A = 
and vanishes for |A| > A* (see the left panel of figure [3]). The latter can be found by 
using the relation p{X) oc linQ{z = A + ie) and the link between Q{z) and B{z). Simple 
reasoning shows that the boundary A* of the domain of existence of eigenvalues is the 
solution of equation B'Jz) = [33] : 



A^ 




-\ — arccoth 

TT 



1 + 



(39) 



This equation simplifies to A* = 2^/^ for /3 <C 1 and to A* = |/3 for /3 3> 1. 

Another important limit of (1371) is that of high density pX^ ^ 1 of points in a large 
box L ^ 1/fco- In this limit the arguments of arctan functions in f l37l) are small and we 
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can put arctanx c^ x. Taking the limit of p\^/ (3' ^ i/k^L ^ we then obtain 

5c(^) = ; + i^l + f^^, pA^»l. (40) 

For pAg below a critical value {pXl)c = 30.3905 the eigenvalue distribution corresponding 
to ( llOl) is asymmetric but bell-shaped, similarly to the case of low density. For 
pXl > (pAo)c, however, the distribution develops a gap: p{X) = for Ai < A < A2, 
where Ai^2 = ^0(^1,2) with Zi^2 being solutions of B'Jz) = (see the right panel of 
figure [3j). In the limit of pAg ^ {p^Dc we have 

A, ~ - ^ - ^ (41) 

'- 27r2 2pA3' ^ ' 

3 ^2/3 ^2 

A. - - 5^(PA2)V3 + —^ + _. (42) 

In figure [4] we compare p{X) following from ( l37l) with the results of numerical 
simulations. We find the Green's function Q^i^) ^^ solving the equation B^[Q^{z)] = z 
numerically and then evaluate the probability distribution of eigenvalues p(A) with the 
help of (IT5I) . When /3 ^ 0, the distribution p{X) tends to the Wigner semi-circle law. 
In contrast, for large /3 > 1 it resembles a Cauchy distribution. A good agreement 
between numerical results and f l37l ) is observed not only for /3 < 1 (similarly to the 
case of sine matrix in section Hj) but for /3 > 1 as well. Note that in contrast to the 
Marchenko-Pastur law f l26l) parameterized by a single parameter /3, the Green's function 
( 1371) and the corresponding probability distribution depend on two parameters /3 and 
pAg. A good agreement between f l37l ) and numerical simulations is obtained at low 
densities pAg < 30 (see figured]). In contrast, at higher densities pAg ^ 30 (not shown) 
the probability distribution following from ( l37j) develops a gap that is not present in 
numerical results. Interestingly, this gap in the probability distribution appears at the 
same density pAg '^ 30 for all /3. 

6. Eigenvalue distribution of cose + i sine matrix 

The matrices C and S can be combined in a single complex non-Hermitian matrix: 
C + \{S — I). The theory of free random variables |^ allows one to study the statistical 
distribution of the complex eigenvalues of this matrix based on the properties of the 
matrices C and S that we considered in the previous sections [3lj. This, however, 
requires asymptotic freeness of C and S. Unfortunately, the matrices S and C defined 
by ( I2T]) and ( l30l) through the same set of points {r^} turn out to be not asymptotically 
free. We therefore start our study of non-Hermitian Euclidean random matrices by the 
case of a matrix X = (7 + 1(5^' — I), where two different and independent sets of points 
{r^} and {r^} are used to define the real and imaginary parts of X\ 
r -(A X xCos(fc^^V2r^ 

f^Ol^i ~ A j I 

_ sin(fco|r:-r;.|) 

'^ ~ kJr' r'\ ' ^ ^ 
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Figure 4. Probability density of eigenvalues of a square N x N Euclidean matrix 
C with elements Cij = (1 — Sij) cos{ko\ri — rj|)//co|r^ — rj|, where the N points 
r^ are randomly chosen inside a 3D cube of side L. Numerical results (blue solid 
lines) obtained for TV = 10^ after averaging over 10 realizations are compared to our 
equation (|37|) (red dashed lines) with P = 2.SN/{koL)'^ for several densities p of points 
(Ao = 27r/A:o). 



The matrix X defined in this way is similar to the matrix G defined in the introduction 
except that it has no correlation between its real and imaginary parts. Using the 
definition of asymptotic freeness [27l [28] it is easy to check that the matrices C and S' are 
asymptotically free, in agreement with the intuitive definition of freeness as statistical 
independence. One can easily show that for the same reason as the one that ensured 
positiveness of the eigenvalues of the matrix S in section [U the complex eigenvalues A 
of the matrix X obey ImA > — 1. 

For non-Hermitian matrices, the Green's function loses its analyticity inside 
two-dimensional domains ('islands') on the complex plane, instead of segments of 
the real axis in the Hermitian case. In [34] Jarosz and Nowak provide a simple 
algorithm, based on the algebra of quaternions, to calculate the non-holomorphic 
Green's function Gxi^) and the correlator of left \Li) and right \Ri) eigenvectors 



[35] Cj^{z) = —{n / N) {J2i=i{Li\Li) {Ri\Ri)S{z — A^)) inside these domains for any non- 
Hermitian matrix of the form X = Hi + iH2^ where Hi and H2 are two asymptotically 
free Hermitian matrices with known Blue functions. In our case. Hi = C and H2 = S^—I. 
In the limit of /3 <C 1, the Blue functions are Bi{z) = (3z + 1/z (section [5]) and 
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B2{z) = l/{l — [3z) — l + l/z (section H]). Qx{z) andC^(2:) can be then found analytically: 

y , 1 



gj^{z = x + iy) = 



X 

2^ 



P{l + y) 2 + y^ 



(44) 



Cjiiz 



x + w)= [y^ 



X 



1 

+ 4 



y 



1 



/3{l + y) 2 + y 



1 



/3(l + y)(2 + y)- ^^^^ 

The correlator (H5l) must vanish on the borderline of the eigenvalue domains. We 
therefore readily obtain an equation for the borderline of the domain of existence of 
eigenvalues of X on the complex plane: 

.-+r^-^V- .. 'l_ . ^0, (46) 



1 + y 2 + y J (l + y){2 + y) 
where x = ReA and y = ImA. The probabihty density inside this domain is 



P(.x,y) = TT [dxRegji{x,y) -dylmg^{x,: 



1 

4n 



1 1 

p^W+y? 



(47) 
If we use this equation 



(2 + |/)2j 

A better model for the Blue function of the matrix C is 
instead of Bi{z) = (3z + 1/z above, analytic calculation becomes impossible but we can 
still compute Gxi^) ^^d C^(^) numerically. The resulting borderline of the eigenvalue 
domain is shown in figure [5] (dashed lines) together with the eigenvalue distribution of 
the matrix X = C + i{S^ — I) found by the numerical diagonalization of a set of 10^ x 10^ 
random matrices. At the smallest density considered pXl = 0.01, the borderline found 
using f l38l ) is very close to f l46l ). At higher densities the former describes numerical 
results much better than (l46l) . 

Equation f l46l ) predicts a splitting of the eigenvalue domain in two parts at /3 = 8. 
The more accurate calculation using ( l38l) makes a similar prediction (see the lower right 
panel of figure Ej). However, the eigenvalues of the matrix X do not show such a splitting 
and form an 'inverted T' distribution on the complex plane instead. This is due to the 
fact that the Marchenko-Pastur law f l26l ) fails to describe the eigenvalue distribution of 
the matrix S'' at /3 > 1 and hence the Blue function 1/(1 — ^z) + 1/z that we assumed 
for S' is not a good approximation anymore. 

It is worthwhile to note that large random non-Hermitian matrices similar to our 
matrix X were considered previously by Haake et al [36j (with the help of the replica 
trick), Lehmann et al [37] (using the supersymmetry method) and Janik et al [29] 
(using the free probability theory). These authors studied matrices of the form /f + icr, 
where H was an Hermitian matrix with random elements obeying Gaussian statistics, 
r was a Wishart random matrix [i.e. a matrix of the form fl25l)]. and c was a real 
number controlling the 'degree of non-Hermiticity' of the matrix. The splitting of the 
domain of existence of eigenvalues in two parts was observed when c was increased. 
This is diflFerent from our matrix X that has elements with equal variances f3/N of real 
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Figure 5. Density plot of the logarithm of the probability density of eigenvalues A 
of a square N x N Euclidean matrix X with elements Xij = (1 — Sij)[cos{ko\ri — 
rj\)/ko\Ti — Tj\ + isin(/co|r^ — r^|)//co|r^ — r^|] at 4 different densities p of points r^, r/ 
per wavelength Aq = 27r//co cube. 2N = 2 x 10^ points r^ and r'- {i = 1, . . . ,7V) are 
randomly chosen inside a 3D cube; the probability distributions are estimated from 10 
realizations of {r^} and {r^}. Dashed lines show the domain of existence of eigenvalues 
following from the free probability theory. 



and imaginary parts (hence always the same degree of non-Hermiticity) but that still 
exhibits the splitting of the eigenvalue domain when /3 is increased. 



7. Eigenvalue distribution of the complex expc matrix 

By analogy with the cardinal sine and cosine functions, a 'cardinal complex exponent' 
function can be defined as f{x) = exp(ix)/x. The Euclidean random matrix G 
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corresponding to this function has elements 

G., = /(r.-r,) = (l-«,)^2^1^^. (48) 

^0|m — Tj\ 

This matrix has a particular importance in the problem of wave scattering by an 
ensemble of A^ point-like scatterers: indeed, as we noted already in the introduction, 
each element of the matrix (7 is a Green's function of the scalar Helmholtz equation ([T]). 
Although the matrix G is similar to the matrix X considered in the previous section, 
the analytic study of its properties is much more involved. On the one hand, similarly 
to the eigenvalues of X, the eigenvalues of G obey ImA > —1. On the other hand, 
correlations that arise between the real and imaginary parts of G due to the presence of 
the same set of points {r^} in both ReG = C and ImG = S' — I, do not permit to take 
full advantage of the free probability approach described in section [6l Another way to 
deal with non-Hermitian matrices is to double the size of the space and to manipulate 
Hermitian matrices of size 2N x 2N (in the 'quaternion' space [29j or in the 'chiral' space 
[38]). Due to technical difficulties, however, this approach can be readily put in practice 
only in certain special cases like, e.g., in the case of circularly invariant distributions 

p{x) = p{\x\) m- 

Despite the differences between the matrices G and X, a comparison of their 
eigenvalue distributions appears to be quite useful. Numerical calculations show that, 
roughly speaking, the eigenvalues of G are concentrated within a circle on the complex 
plane (see figure [6l figure [71 and figure [8]). The same circular shape of the domain 
of existence of eigenvalues is characteristic for the matrix X in the limit of /3 -^ 0. 
The radius of the circle is ^/2^ and the position of its center is Xq = 0, ?/o = |/3 as 
can be seen from ( l46l) by assuming y <C 1 in the denominators. To derive this result 
in a more rigorous way, we substitute the equation of a circle, x = ry/2^co8(j) and 
y = rv^2^sin(/) + |/i/3, in ( l46l) . The l.h.s. of the resulting equation is then expanded 
in orders of /3 and coefficients Cjy in front of consecutive powers of /3 are analyzed. The 
coefficient Cq in front of /3^ is zero whatever r and h. The coefficient in front of /3^ is 
Ci = 2(r^ — 1). Equation ( l46l) is therefore obeyed up to the linear order in /3 for r = 1 
and any h. The next term is of the order /3^/^ and it cannot be put to zero simply 
by adjusting h because C3/2 depends on 0: Cs/2{(t>) = A/2sin0(/i + 2 cos 20). We thus 
search for h that minimizes the integral J^^ d(j) Cy2{^) = 27r(/i^ — 2h + 2). This yields 
h = 1. The density of eigenvalues inside the circle is not homogeneous: expanding 
( l47ll around y = yo and taking the limit /3 ^ we obtain p{x^y) 0^ {1 — y)/2Ti(5. 
In figure [6] we superimpose the circle of radius ^/2^ centered at Xq = 0, ?/o = |/3 on 
the eigenvalue distribution of G for small /3. The circle describes the boundary of the 
eigenvalue distribution remarkably well. We thus conclude that in the limit of /3 <C 1, 
the domains of existence of eigenvalues of the matrices G and X are very similar. In 
addition to the eigenvalues inside the circle, G has eigenvalues that follow the spirals 
corresponding to the eigenvalues G12 and — G12 of a 2 x 2 matrix G. Interestingly, the 
spirals are quite robust and survive at all densities (see figure [6l figure [71 and figure [8]) . 

Because the matrices S and C studied in previous sections represent the imaginary 
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Figure 6. Left column: density plot of the logarithm of the probability density 
of eigenvalues A of a square N x N Euclidean matrix G with elements Gij = 
(1 — Sij) ex.p{iko\Ti — Tj\)/ko\Ti — Tj\ at low densities pAg = 0.01 (first row) and 0.1 
(second row), Aq = 27r/A:o. A^ = 10^ points r^ are randomly chosen inside a 3D cube. 
Dashed circles are centered at (0, ^/3) and have radii y/2p. Eigenvalues of a 2 x 2 
matrix would lie on dashed spirals. Central column: marginal probability density of 
the real part of A compared to our equation (|37|) with /3 replaced by ^/3 (dashed red 
line). Right column: marginal probability density of the imaginary part of A compared 
to the Marchenko-Pastur law (|26]) with A replaced by ImA + 1 and (3 replaced by ^/3 
(dashed red line). 



and real parts of the matrix G, respectively, one might expect some links between 
the probability distributions of eigenvalues of S and C and the marginal probability 
distributions of the real and imaginary parts of the eigenvalues of G. And indeed, 
we see from the central and right columns of figure [6] that the marginal probability 
distributions p{ReX) and p{linX) are nicely described by f l37l) and f l26ll . respectively, 
with /3 replaced by |/3. This suggests an interesting interpretation of the Marchenko- 
Pastur law f l26l ): it can be seen as a projection of a two-dimensional distribution p{x^ y) 
of complex eigenvalues x + i^ on the imaginary axis y, provided that p{x^ y) is different 
from zero only inside a circle of radius 2^/^ centered at (0,/3) and that p{x^y) oc 1/y 
inside the circle. p{x^ y) being independent of x and decaying monotonically with y is 
consistent with the result that we obtained for the matrix X using the free random 
variable theory in the limit of /3, pAg -^ 0. 

When we increase (5 but keep the density relatively low (pAg < 30, see below), the 
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Figure 7. Same as figure [6] but for intermediate densities pAg = 1 (first row) and 10 
(second row). Dastied circles are centered at (0, ^f3) and have radii R given by (|49|) . 



cloud of eigenvalues grows but keeps its circular shape (see figure [7j) □• The distribution 
of eigenvalues acquires an important asymmetry: the eigenvalues are 'attracted' by the 
axis ImA = — 1. Interestingly, whereas the Marchenko-Pastur law ceases to describe the 
marginal distribution of ImA when |/3 becomes larger than unity, the region of validity 
of f l37l ) for the distribution of ReA is wider: as we show in figure [71 f l37j) continues to yield 
reasonable results even for |/3 > 1. As can be seen from figure [71 even at |/3 > 1 the 
borderline of the eigenvalues' domain is still roughly a circle. More accurate inspection 
reveals that this circle is still centered at (0, |/3) even for /3 ;» 1. It touches the line 
ImA = — 1 that it cannot cross. Its radius is, therefore, roughly |/3 and not ^/2^ as in 
the limit of small /3. To extrapolate between the limits of small and large /3 we propose 
the following empirical expression for the radius R of the eigenvalue domain: 

i?2^2/5+(^0 . (49) 

For /3 <C 1, the second term of this equation is negligible and we recover R = \/2^. For 
the parameters of figure [6l for example, a circle of radius R given by f l49l) is virtually 
indistinguishable from the circle of radius ^/2^ shown in the figure. At larger /3 the 

+ Because we present results at a fixed N = 10^, increasing /3 is achieved by increasing the density 
pXq. However, by repeating the analysis at A^ = 10^ and A^ = 5 x 10^ we checked that the distributions 
presented in figure [71 change only slightly when N and /)Aq are varied to keep /3 constant. 
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Figure 8. Same as figure [6] and figure [71 but for high densities pAg = 50 (left column) 
and 100 (right column); dashed circles as in figure [71 Note a hole that develops on the 
left from ReA = near the real axis and the corresponding gap in the analytic result 
for the marginal distributions of ReA (dashed red lines). Marginal distributions of ImA 
are not shown. 



second term in f l49l) starts to play a role and dominates for /3 :::> 1. As we show in 
figure U\ ( l49l) gives a good idea of the part of the complex plane where the eigenvalues 
of the matrix G are concentrated. 

At high densities pAg ^ 30, a 'hole' appears in the eigenvalue distribution that 
otherwise still preserves its overall circular structure (see figure [8j) [j. Interesting enough, 
this hole is not accompanied by any visible signatures in the marginal distributions 
p(ReA) and p{linX). However, the analytic result f l37l) develops a gap precisely at 
the same density pAg ^ 30 and at the same position at which the hole appears on 
the complex plane. This suggests that even though f l37l) does not provide a correct 

* By repeating calculations with N = 10^ and TV = 5 x 10^ we found that the density pAg at which 
the hole appears in the eigenvalue distribution is roughly independent of f3. 
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Figure 9. Marginal distribution of the imaginary part of the eigenvalues of the matrix 
G computed numerically at densities pX^ = 1, 10, 20, 40, 60 and 100 (curves from top 
to bottom) for N = 10^ are compared with the asymptotic law l/F shown by the 
dashed hue. 



description of the marginal distribution p{ReX) at such high densities, it still reflects 
some relevant properties of the distribution of complex eigenvalues A. Note that at high 
densities pXl the eigenvalue distribution is concentrated near the axis ImA = — 1 and the 
parts of the distribution corresponding to ImA ::^ 1 in figure [8] are visible only thanks 
to the logarithmic scale of the plot. 

Finally, we study the marginal distribution of ImA. It has been given special 
attention previously because, under certain assumptions, it was shown to give the 
distribution of 'decay rates' F = ImA + 1 of quasi-modes in an open random medium 
[inilll]. When |/3 > 1, p{linX) does not follow the Marchenko-Pastur law anymore (see 
figure[7j). Based on the results of numerical simulations, Pinheiro et al. [11] claimed that 
at high densities p\^ the marginal distribution pilmX) exhibits a universal l/F decay. 
Our analysis summarized in figure [9] confirms that such a decay is present, even though 
it seems to speed up slightly when the density is increased. In certain applications 
of random matrix theory to wave propagation in random media and, in particular, 
in problems related to Anderson localization (see section 18.21) and random lasing (see 
section 18.31) , a special role is played by the eigenvalue of G that have the smallest or 
the largest imaginary part. Both min(ImA) and max(ImA) are random variables. Let 
us first consider min(ImA). As can be seen from figure [6] and figure [71 at moderate 
densities pAg ^ 10, min(ImA) is due to the lower spiral emerging from the 'bulk' of 
the distribution. Eigenfunctions of G corresponding to spirals are localized on pairs 
of nearby points and the eigenvalues can be found by considering a 2 x 2 matrix G. 
For two points at a distance Ar we find the eigenvalues Ai^2 = ±exp(i/coAr)//coAr, 
with A2 corresponding to the lower spiral. The smallest values of ImA are achieved 
for small distances Ar when we can approximately write ImA2 = — sin(A;oAr)/A;oAr c::^ 
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— 1 + {koAry/6. Hence, the statistical distribution of min(ImA) is directly related 
to the statistical distribution p(Armin) of the minimal distance Armin between any 2 
points among A^ points in the volume V. The distribution p(Armin) can be constructed 
as follows. Let us choose an arbitrary point i. The probability that another point j is 
located in a spherical shell of radius r and thickness dr around the first is pi = Aur^dr /V . 
For r to be the minimal distance Armin we have to require that all other N — 2 points 
are outside the sphere of radius r [probability p2 = {I — 4:7Tr^/3V)^~'^] and that the 
distances between the remaining (A^ ~ 1)^ pairs of points not including the point i 
exceed r [probability ps = (1 — 47rr^/31/)^^~^^ ]. The probability that r is the minimum 
distance between any 2 points is then equal to the number of possibilities A^(A^ — 1) to 
choose the two points i and j, times Pi x P2 x Ps- The probability density is then 

piAr^in) = N{N - 1) (^^^^ [l - ^^^ • (50) 

This distribution is normalized to 1 if we assume that the volume V is spherical (radius 
Rq) and that Armm can vary from to i?o. Because min(ImA) = — 1 + fcoAr^-^^/G, its 
probability density is equal to p[Ar^i^ = ^6(min(ImA) + l)/fco] x [dArmin/dmin(ImA)]. 
In particular, the first moment of this distribution in the limit of A^ ^ oc is 

(mindmA)) = -1 + '-^^^^ x (,,3,V)V3 ' <"> 

We compare this result with numerical simulations in figure [10] (left panel) and find 
good agreement for densities pX^ < 10. At higher densities, (min(ImA)) is smaller 
than predicted by (I5TI) . signaling that min(ImA) is not dominated by the eigenvalues 
corresponding to eigenfunctions localized on pairs of points anymore. 

Similarly to min(ImA), max(ImA) is dominated by the second spiral branch of the 
eigenvalue distribution for /3 < 0.3 (see figure [6]). At larger /3, max(ImA) belongs to 
the bulk of the eigenvalue distribution (see figure [7] and figure [8]) . As follows from our 
analysis, the distribution of complex eigenvalues A of the matrix G occupies a circular 
domain of radius R given by f l49l) . centered at |/3. It follows then that 

(max(ImA)) ^ ^f3 + R. (52) 

And indeed, this approximate expression describes numerical results quite reasonably 
(see the solid line in the right panel of figure fTOl) . even though a closer inspection reveals 
that it overestimates (max(ImA)) at large /3. Further work is needed to find a more 
accurate expression for (max(ImA)). 

8. Applications 

We have already mentioned in the introduction that the Euclidean random matrices S, 
C and G studied in this paper are encountered in several physical problems. In this 
section we briefly discuss a number of such problems and show how our results can help 
to advance their understanding. 
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Figure 10. Mean minimum (left) and maximum (right) values of the imaginary 
part of eigenvalues A of the matrix G for three different matrix sizes N (symbols). 
(Min(ImA)) + 1 is approximately 2.3[N{pXq)]~'^^^ [solid line in the left panel and (|5T]) ] 
for pAg ^10 and decays faster at higher densities. (Max(ImA)) scales with (3. Different 
values of /3 are obtained by changing the density pAg from 0.01 to 100. The solid line 
in the right panel shows (max(ImA)) = ^/3 -\- R with R given by (|49|) . 



8.1. Cooperative emission of large atomic clouds 

An interesting problem of modern quantum optics is the one in which a single photon 
is stored in a cloud of (cold) atoms. One studies the properties (frequency, direction of 
propagation, etc.) of the photon re-emitted by the cloud at a later time [HI [151 flGlfTT] . 
For A^ two-level atoms (excited state a, ground state b) located at random points r^, 
i = 1, . . . , A, the state of the system at a time t can be written as [ITj 



Mt) 



N 

E 






f3j{t)\hb2 



•M|o) + J]7k(t)|&i&2---M|ik) 



Mill 



(53) 



m<n k 

Here the first sum corresponds to the superposition of states with one atom (atom j) 
in the excited state, all other atoms in the ground state, and zero photons. The second 
sum corresponds to the states in which all atoms are in the ground state, while there is 
a photon in the mode k. Finally, the last sum describes states with atoms m and n in 
the excited state and one virtual photon with 'negative' energy. 

The evolution equation for the vector /3{t) = {/3j{t)} reads [TGlfTTj : 

(3{t) = -Tom + ^^oGf3{t), (54) 

where Fq is the spontaneous decay rate of a single atom and the matrix G is defined by 
( l48l ). According to this equation, a system prepared in the eigenstate described by a 
vector /3{0) decays with a rate Fo(l + ImA) and experiences a frequency shift — FgReA, 
where A is an eigenvalue of the matrix G. Both the decay rate and the frequency shift 
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were studied in [161 [El i^i the limit of a very dense atomic cloud (pAg -^ oc), when the 
summation [G/3{t)]j = X]m=i Gim(^m{t) can be replaced by integration in the last term 
on the r.h.s. of f l54l ). The authors also discussed a useful approximation in which the 
real part of the matrix G is neglected and G is replaced by iS in f l54l ) . 

Although the results of [I6l [IT] are very interesting, atomic clouds of moderate 
density pAg % 1 are readily created in modern laboratories (see, e.g., [22l [2^). It is 
therefore important to extend the analysis of [Ml [IT] to such dilute atomic clouds. The 
present work provides, in fact, such an extension: the distribution of dimensionless decay 
rates V = 1 + ImA is given by the Marchenko-Pastur law f l26l ) with /3 replaced by |/3 and 
the distribution of dimensionless frequency shifts Q = —Re A follows from the analysis 
of section [5] (see also figure [6l figure [71 figure [81). It is important to realize that replacing 
summation by integration in the last term on the r.h.s. of f [54l ) performed in [T6l \T7\ 
is equivalent to averaging this equation over all possible configurations {r^} of atoms. 
It leads, therefore, to the neglect of the statistical nature of the initial problem. In 
contrast, our treatment does not rely on such an averaging and fully accounts for large 
fiuctuations of eigenvalues, typical for situations when light is scattered in a strongly 
disordered environment. As a consequence, the authors of [161 [El find deterministic 
eigenvalues A^, whereas we work with the probability distribution p{X). Our results are 
consistent with those of [lEl [El in the limit of pAg -^ oc and provide a generalization 
of some of them. For example, the authors of [Ml |E1 predict that for pX^ -^ oc, the 
fastest decay rate Fmax is of the order of /3. Our study suggests that the dependence 
of Tjnax on (3 (and not on the density pX^) is a general property valid at any density 
(see the right panel of figure [IHj) as well as it yields a more precise relation between 
Tmax = 1 + max(ImA) and /3 [see f[52l)]. 



8.2. Anderson localization in an open medium 

The phenomenon of Anderson localization is common for all waves in random media 
[24l [25l [26] . It consists in a transition from extended (over the whole available sample 
volume) to exponentially localized eigenstates of a wave (or Schrodinger) equation 
with a randomly fiuctuating dielectric constant (or potential), at a sufficiently strong 
randomness. A paradigm system in which Anderson localization can be studied for 
classical waves is a random arrangement of A^ identical point-like scatterers in a volume 
V. In such an open system of finite size the wave energy can leak to the outside and 
one expects Anderson localization to have an impact on decay of physical observables 
(such as, e.g., the intensity of the wave emerging from the random system). Given a 
simple model for scatterers, the relevant decay rates are related to the imaginary part 
of the eigenvalues A of the non-Hermitian matrix G [TOj . 

Several authors studied the distribution of dimensionless decay rates F = ImA + 1 
in open random media and, in particular, promoted the idea of using its probability 
distribution p{r) as a criterion for Anderson localization [HI [401. More precisely, p{r) 
is expected to decay as 1/F in the localized regime. Our numerical results also exhibit 
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such a behavior (see figure [91), but we cannot claim any relation between it and Anderson 
localization. Indeed, a careful inspection of our results shows that p(r) starts to exhibit 
1/r behavior right after the criterion |/3 < 1 breaks down. On the one hand, for 
resonant point-like scatterers the mean free path can be estimated in the independent 
scattering approximation as ^ = I /pa = kl/Aup with the resonant scattering cross- 
section a = 47r/A;Q. The criterion |/3 = 1 corresponds then to a condition for the optical 
thickness L/l o:^ 9. On the other hand, Anderson localization is expected to take place 
for koi c^ 1 (loffe-Regel criterion [25j) which can be rewritten as pAg c^ 20. We thus see 
that the condition required to observe l/F decay of p{r) {L/£ > 9) does not seem to 
agree with the one expected for the Anderson localization (pAg ^ 20). 

The results that we obtained in the present paper suggest another way of using 
statistics of eigenvalues of G to look at the transition from weak to strong scattering 
and eventually to Anderson localization. First, instead of studying the imaginary part of 
A one can study its real part. At low density pAg <C 1 the distribution ]9(ReA) exhibits a 
transition from the Wigner semi-circle law for /3 ~ L/i <C 1 (see figure [6]) to the Cauchy 
distribution for /3 ^ L/i ^ 1 (see figure [7]) . This transition can be seen as a signature 
of the change of regime of wave scattering from single (for L/i <C 1) to multiple (for 
L/i ^ 1) scattering. Second, an important modification oi p{X) that takes place when 
the density of scatterers is increased is the appearance of a hole in the distribution that 
otherwise occupies a circular domain on the complex plane. The condition pAg c^ 30 
for the appearance of the hole is remarkably close to the condition pAg — 20 expected 
for the Anderson localization transition in the independent scattering approximation. 
A highly speculative conjecture might be that a link exists between the hole in p{X) 
on the complex plane and Anderson localization of waves in an ensemble of point-like 
scatterers. Further work is required to prove or to refute this conjecture. 

8.3. Random lasers and optical instabilities 

The eigenvalues of the matrix G that have the smallest imaginary part play a particularly 
important role for understanding of very interesting optical systems called 'random 
lasers'. Random laser is a laser that have no external cavity and in which the feedback 
is provided by the multiple scattering of light J4T[ [42l 143] . One of the minimal models 
to study random lasing is an ensemble of point-like scatterers ('atoms') randomly 
distributed in a volume V = L^ filled with some continuous amplifying medium that 
provides a constant amplification rate Fampi- Lasing starts when Fampi becomes larger 
than the minimum loss rate Fmin = 1 + min(ImA). Therefore, the average value 
of min(ImA) defines the average random laser threshold: (F^^pJ = 1 + (min(ImA)). 
Pinheiro and Sampaio [I9] studied this latter quantity numerically and found a scaling 
law 

1 + (min(ImA)) oc j^2/3^pXiy/3 - (^5) 

They provided a simple interpretation of this result in terms of the diflFusion theory of 
light scattering. 
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The eigenvalues that have the smaUest imaginary part also define the threshold for 
dynamic instabilities in nonlinear random media. In particular, a random arrangement 
of point-like nonlinear scatterers with an intensity-dependent scattering matrix t{I) = 
(— 27ri/fco)[exp(ia/) + 1] was considered by Gremaud and Wellens [SI]. Here / is the 
intensity of light on the scatterer. It was shown that stationary, time-independent 
solutions lose their stability and the system starts to exhibit complex, spontaneous 
dynamic behavior when the nonlinear coefficient a exceeds a critical value ainst- 
The average value of the instability threshold was found to scale as (o^inst) oc [1 + 
(min(ImA))]^/^, with 

This result can be explained by considering the spiral branches of the statistical 
distribution of A on the complex plane (see the dashed spirals in figure [6l figure [7] 
and figure [8j) ; it is not related to the diffusion of light in the bulk of the random sample 
but originates from sub-radiant states localized on pairs of mutually close scatterers 

As follows from the aforesaid, the results f l55l) and f l56l ) that are supposed to 
coincide, not only differ by a factor (pAg)"^/^ but they are given different physical 
interpretations as well. The analysis that we performed in this work allows us to resolve 
this controversy and to identify the result f l56l) of Gremaud and Wellens [21] as the 
correct one. Moreover, not only we are able to derive an analytical expression f lSTI) 
that agrees with f l56l) and contains the precise numerical coefficient, but also the full 
distribution function of min(ImA) — and hence of the instability threshold o^inst — 
follows from our result f l50l) . 

A random laser different from that considered in [19] is the one in which the 
amplification is provided by the point scatterers themselves and not by the medium 
in between them. Amplification and scattering in such a system are not independent 
anymore and cannot be tuned at will. The eigenvalues A of G governing the laser 
threshold will now depend on the specific amplification scheme. Curiously, for the 
simplest physical pumping mechanism we can think of (incoherent pump), the laser 
threshold will be determined by the eigenvalues having the largest imaginary part pO] . 
This provides a direct physical application for the results that we show in the right panel 
of figure [TOl and will be discussed in detail elsewhere [20] . 

9. Conclusion 

In this work we studied eigenvalue distributions of certain Euclidean random matrices 
that appear in the context of wave propagation in random media. In particular, 
we considered large N x N real symmetric matrices S and C with elements Sij = 
sin(/co|r^ — rj|)//co|r^ — rj| and Cij = (1 — Sij) cos(/co|r^ — rj|)//co|r^ — rj|, respectively, as 
well as the non-Hermitian matrix G = C + i{S — I).N points r^ were chosen randomly 
in a three-dimensional cube of side L with density p = N/L^. For the three random 
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matrices under study, the two important parameters of the eigenvalue distributions p{X) 
are /3 = 2.8N/{koLy and the number of points per wavelength cube pAg. /3 is equal to 
the variance of eigenvalues A of both S and C in the limit of koL -^ oc. 

In the low-density limit pAg <C 1 and for /3 < 1, the distributions of eigenvalues 
of Hermitian matrices S and C are parameterized uniquely by /3: the distribution of 
eigenvalues of S is given by the Marchenko-Pastur law f l26l ). whereas the distribution 
of eigenvalues of C can be deduced from the Blue function f l38l ) that we derived in this 
paper. For /3 > 1 the Marchenko-Pastur does not apply to S anymore, but out equation 
( l38l ) still works for C as long as pAg is small enough. As /3 increases, f l38l ) describes a 
transition from the Wigner semi-circle law (at /3 <C 1) to the Cauchy distribution (at 
/3 ;» 1). At high densities pAg > 1 the more complete expression ( 1371) that we derived 
for the Blue function of the matrix C applies. It is in good agreement with numerical 
simulations until pAg ^ 30 where it predicts the appearance of a gap in p(A), which is 
not observed in the numerical data. 

The eigenvalue distribution of the non-Hermitian matrix G has a circular structure 
on the complex plane. At /3 <C 1, the eigenvalues are confined to a circle of radius \/2^ 
centered at (0, |/3). At larger /3, the distribution becomes strongly asymmetric, with 
much stronger weight of eigenvalues with imaginary parts close to —1. Our numerical 
results show that the domain of existence of eigenvalues is still approximately a circle 
centered at (0, |/3). We proposed an empirical expression for its radius i?^ ^ 2^ + (|^)^. 
At high densities pAg > 30, a hole appears in the distribution p{X) on the complex 
plane. The density at which the hole appears seems to be roughly independent of (3. 
The marginal probability distribution of Re A is described by our equation f l37l) at all /3, 
provided that pAg < 30. The marginal distribution of ImA follows the Marchenko-Pastur 
law f l26l ) for |/3 < 1 and decays as l/(ImA + 1) at larger /3. 

Finally, we studied a model matrix X = C + i{S' — I) in which two independent 
ensembles of points {r^} and {r^} were used to generate matrices C and S\ The matrices 
C and S' are asymptotically free and the distribution of eigenvalues of X at /3 < 1 can 
be found using the approach developed by Jarosz and Nowak [3lj based on the theory 
of free random variables. The distribution of eigenvalues shows an interesting transition 
from a circular shape at /3 <C 1 to a triangular shape at /3 ^ 1, and then to an 'inverted 
T' shape for /3 > 1. 
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